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PROCESSING OR COMPRESSING N-DIMENSIONAL SIGNALS WITH WARPED WAVELET PACKETS 
AND BADLBTS 



BACKGROUND OF INVENTION 

5 The invention relates generally to a n-dimensionai signal processing method, apparatus and computer pro- 
gram, and in particular to a method, apparatus and computer program useful for processing n-dimensional 
signals, such as one-dimensional signals, two-dimensional images or three-dimensional data or video image 



The invention is particularly pertinent to the field of signal data processing and compression. Signal 
10 compression is a process which encodes signals for storage or transmission over a communication channel, 
with fewer bits than what is used by the uncoded signal The goal is to reduce the amount of degradation 
introduced by such an encoding, for a given data rate. The invention is also relevant to signal restoration or 
feature extraction for pattern recognition. 

In signal processing, efficient procedures often require to compute a stable signal representation which 
15 provides precise signal approximations with few non-zero coefficients. Signal compression applications 
are then implemented with quantization and entropy coding procedures. At high compression rates, it has 
been shown in S. Mallat and F. Falzon, "Analysis of low bit rate image transform coding," IEEE Trans. 
Signal Processing, voL 46, pp. 1027-1042, 1998, that die efficiency of a compression algorithm essentially 
depends upon the ability to construct a precise signal approximation from few non-zero coefficients in the 
20 representation. 

The stability requirement of the signal representation has motivated the use of bases and in particular 
orthogonal bases. The signal is then represented by its inner products with the different vectors of the basis. 
A sparse representation is obtained by setting to zero the coefficients of smallest amplitude. During die last 
twenty yeare, different signal representations have been constructed, with fast procedures which decompose 

25 the signal in a separable basis. Block transforms and in particular block cosine bases have found important 
applications in signal and image processing. The JPEG still image coding standard is an application which 
quantizes and Huffman encodes the block cosine coefficients of an image. More recently, separable wavelet 
bases which compute local image variations at different scales, have been shown to provide a sparser image 
representation, which therefore improves the applications. These bases are particular instances of wavelet 

30 packet bases, in R. Coifinan, Y. Meyer, M. Wickerhauser, 'Method and apparatus for encoding and decoding 
using wavelet-packets", US Patent 5,526,299. The JPEG image compression standard has been replaced by 
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the JPEG-2000 standard which quantizes and encodes the image coefficients in a separable wavelet basis: 
"JPEG 2000, ISO/EEC 15444-1:2000 " 2000. Wavelet and wavelet packet bases are also used to compress 
one-dimensional signals, including medical signals such as electro-cardiogram (ECG) recordings, as in M. 
Hilton, J. Xu, Z. Xiong, "Wavelet and wavelet packet compression of electrocardiograms", IEEE Trans. 
5 Biomed Eng., voL 44, pp. 394- 402, May 1997. Decomposition in three dimensional wavelet bases 
are also used in video image compression, in S. Li and Y-Q. Zhang, in 'Three-Dimensional Embedded 
Subband Coding with Optimized Truncation (3-D ESCOT)", Applied and Computational Harmonic Anal- 
ysis 10, 290-315 (2001), where a video sequence is decomposed with three dimensional wavelet transform 
performed along motion threads in time. 
10 Signal restoration of sparse signal representations has been developed by thresholding the wavelet co- 
efficients of noisy signals in D. Donoho and I. Johnstone, Ideal spatial adaptation via wavelet shrinkage" 
Biometrika, vol. 81, pp. 425-455, December 1994. Applications of wavelet packet bases to deconvolution 
of signals are also presented in J. Kalifa, S. Mallat, "Minimax restoration and deconvolution", in Bayesian 
inference in wavelet based models, ed. P. Muller and B. Vidakovic, Springer-Verlag, 1999, Constructing 
15 sparse representations is also important to extract features for pattern recognition. This has important ap- 
plications to content based signal indexing and retrieval from digital multimedia libraries and databases. 
Feature vectors using histograms of wavelet coefficients are used in M. K. Mandal and T. Aboumasr, "Fast 
wavelet histogram techniques for image indexing", Computer Vision and Image Understanding, vol 75, no. 
1/2, pp. 99-110, August 1999. 
20 The main limitation of bases such as block cosine bases, wavelet bases or more generally wavelet packet 
bases, currently used for signal representation, is that these bases are composed of vectors having a fixed 
geometry which is not adapted to the geometry of signal structures. For one-dimensional signals such as 
ECG, which are quasi-periodic, adapting the basis to the varying period allows one to take advantage of the 
redundancy due to tie existence of a periodicity in the signal In images, edges often correspond to piece- 
25 wise regular curves which are therefore geometrically regular. In higher dimensional signals, such as video 
sequences, edges and singularities often belong to manifolds that are also geometrically regular. Construct- 
ing bases that take advantage of this geometrical regularity can considerably improve the efficiency of signal 
representations and hence improve applications such as compression, restoration and feature extraction. 
In E. Le Pennec and S. Maflat, 4t Method and apparatus for processing or compressing n-dimensional 
30 signals by foveal filtering along trajectories", US Patent application Serial No. 09/834,587, filed April 13, 
2001, and in E. Le Pennec, S. Mallat, "Image Compression with Geometrical Wavelets", Proceedings of 
International Cpn£ on Image Processing, Vancouver, September 2000, part of the signal information is rep- 
resented with wavelet foveal filters that follow foveal trajectories adapted to the geometry of the signal. The 
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wavelet fovea! coefficients are then decorrelated with linear operators that compute bandelet coefficients. 
The edgeprint representation of Dragottia and Vetterli, in "Footprints and edgeprints for image denoising 
and compression", Proceedings of the International Conference on Image Processing, Thessaloniki, Octo- 
ber 2001, use a similar strategy with footprint wavelet vectors that follow edges computed from the image. 
5 Fovea! bandelets and edgeprints do not provide a complete signal representation, and it is therefore neces- 
sary to incorporate a residua! signal to reconstruct the original signal, which is a source of inefficiency for 
data compression and restoration applications. > 

Accordingly, there exists a need in the art for improving signal processing, by computing sparse rep- 
resentations by taking advantage of the signal geometrical regularity, from which one can reconstruct pre- 
10 cise signal approximations with fast and numerically stable procedures and apply it to signal compression, 
restoration and pattern recognition. 

SUMMARY OF THE INVENTION 

It is a primary 'object of this invention to devise a method and means to construct a sparse and stable warped 
wavelet packet representation of n-dimensional signals by taking advantage of the regularity of their geo- 

15 metrical structures. A further object is to a compute bandelet representation from the warped wavelet packet 
representation,- with an efficient bandeletisation adapted to the signal geometry. It is yet another object of 
this invention to build a system that compresses signals by quantizing and encoding the coefficients of this 
sparse bandelet representation. Another object of this invention is to restore signals by processing the coef- 
ficients of this bandelet representation. Another object of this invention is to use the bandelet representation 

20 for signal feature extraction for pattern recognition systems. 

The invention includes a warped wavelet packet processor that computes an n-dimensional warped 
wavelet packet transform including warped wavelet packet coefficients and wavelet packet warping grids, 
from an n-dimensional digital input signal, wherein n is a positive integer. It comprises the steps of provid- 
ing an n-dimensional warped signal including n-dimensional warped coefficients and n-dimensional signal 

25 warping grids; and computing said warped wavelet packet transform of said warped signal, with a binary 
tree where each node performs a one-dimensional warped subband processing along a particular dimension 
<i, with 1 < d < n. In dimension n > 2, said warped subband processing have a phase alignment coherent 
with said signal warping grids. 

The invention also includes an inverse warped wavelet packet processor that computes an n-dimensional 

30 digital output signal from an n-dimensional warped wavelet packet transform* It comprises the steps 
of: computing a warped signal including n-dimensional warped coefficients and n-dimensional signal 
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warping grids, from said n-diroensional warped wavelet packet transform, with a binary tree where each 
node performs a one-dimensional inverse warped subband processing along a particular dimension d, with 
1 <d <n; and computing said digital output signal from said warped signal with an inverse signal warping. 
The sampling grid of the output signal is identical to the sampling grid of the input signal. 
5 As opposed to prior art wavelet packet processors, as in R. Coifman, Y. Meyer, M Wickerhauser, 
"Method and apparatus for encoding and decoding using wavelet-packets" US Patent 5,52639, or to J. Li, 
S JM. Lei, "Arbitrary shape wavelet transform with phase alignment", US Patent 6,233,357, or to A, Merlins, 
Image compression via edge-based wavelet transform," Opt Eng., vol. 38, no. 6, pp. 991-1000, 1999, 
the subband filtering is not performed on the input signal sampling grid, but on different sampling grids, 
10 . called warping grids, that are typically adapted to the geometrical signal properties in different regions. As 
opposed to the subband filtering used in the wavelet packet procedures referenced above, a warped subband 
filtering is not computed with a convolution operator but with space varying filters, in order to handle the 
nonuniform structure of the warping grids. For three-dimensional video image sequences, as opposed to 3D- 
ESCOT method in J. Xu, Z. Xiong, S. Li and Y-Q. Zhang, < Three-Dimensional Embedded Subband Coding 
15 with Optimized Truncation (3-D ESCOT)", Applied and Computational-Harmonic Analysis 10, 290-315 
(2001), the warping grid regions are not reduced to motion threads in time, and the warped subband pro- 
cessing satisfies a phase alignment constraint coherent with the warping grids, which is not satisfied by the 
3D-ESCOT method. With this phase alignment property, warped wavelet packet coefficients have the same 
geometrical regularity as the input signal This is particularly interesting when a bandeletisation module is 
20 located downstream of the warped wavelet packet processor. 

To adapt the signal representation to the geometrical signal structures, the invention comprises a geo- 
metrical sampling module that computes said signal warping grids from a set of parameters called a warping 
geometry. The warping grids are typically computed to follow the directions in which the signal has regular 
geometrical variations. In an exemplary implementation, the warping geometry includes region parameters 
25 that specify a partition of the signal support into several regions, and includes deformation parameters that 
define a geometrical deformation function which specify the position of sampling points in each of said 
region. The signal support is divided into regions in which the signal has typically uniform geometrical 
properties so that one can adapt the warping grid to all structures in the region. For signals that are nearly 
periodic, the warping grid can adapt the sampling interval to the varying period, to obtain a nearly exactly 
30 periodic signal, which is taken advantage of, by the subsequent bandeletisation module. 

Warped wavelet packet coefficients are computed with warping grids that typically follow directions in 
which the signals has regular variations. Because of the phase alignment property of warped subband pro- 
cessors, warped wavelet packet coefficients inherit the regularity of the signal in these directions. The raven- 
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tion preferably includes a bandeletisarion module that yields a sparse representation by decorrelating said 
warped wavelet packet coefficients, by applying invertible one-dimensional linear operators along selected 
directions of said warped wavelet packet coefficients. For regions in which the signal is nearly periodic 
along particular directions, the bandeletisarion module performs a periodic decorrelation, that takes advan- 
5 tage of the redundancy introduced by the existence of a quasi-periodicity. The resulting bandelet coefficients 
are decomposition coefficients in a basis composed of warped bandelet vectors. The one-dimensional linear 
decorrelation operators, can be chosen to be a cosine transform or a one-dimensional wavelet packet trans- 
form or a warped wavelet packet transform. The n-dimensional signal is then represented by its bandelet 
coefficients and the parameters of the warping geometry that specify the warping grids in each signal region, 
10 The invention also includes an inverse bandeletisarion module that computes an n-dimensional warped 
wavelet packet transform from a warping geometry and bandelet coefficients. It comprises the steps of 
computing wavelet packet warping grids from said warping geometry, and computing warped wavelet packet 
coefficients by applying inverse one-dimensional linear operators along selected directions on said bandelet 
coefficients. 

15 ftTopposedto'the method in'E.Le Pennec and S. Mallat, in "Method and apparatus for processing or 

compressing n-dimensional signals by fbveal filtering along trajectories**, US Patent application Serial No. 
09/834,587, filed April 13, 2001, and in E. Le Pennec, S. Mallat, "Image Compression with Geometrical 
Wavelets", Proceedings of International Conf. on Image Processing, Vancouver, September 2000, bandelet 
coefficients are not computed from fovea] coefficients but from warped wavelet packet coefficients. As a 

20 consequence, no residual signal is needed to reconstruct a precise signal approximation from bandelets of 
warped wavelets as opposed to the above referenced method. 

Signal processing procedures are efficiently implemented in a warped wavelet packet bandelet represen- 
tation because of the ability to provide sparse and accurate representations by setting their smallest coeffi- 
cients to zero. The invention comprises a processor compressing n-dimensional signals that quantizes the 

25 bandelet coefficients and encodes these quantized bandelet coefficients with the warping geometry to obtain 
a multiplexed data stream suitable for storage in a storage medium or for transmission over a transmission 
medium. The invention also comprises a processor that restores an input signal by applying a restoration 
processor on the bandelet coefficients and the warping- geometry to compute an output signal from these 
restored coefficients. The invention also comprises a processor that computes a signal feature vector from 

30 the warping geometry and bandelet coefficients, for partem recognition applications including content based 
signal indexing or retrieval and signal matching. 
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BRIEF DESCRIPTION OF THE DRAWINGS 

The foregoing and other objects of this invention, the various features thereof as well as the invention itself 
may be more fully understood from the following description, when read together with the accompanying 
drawings in which: 

5 FIG. 1 shows, in block diagram form, an exemplary embodiment of the invention which computes a 
warped wavelet packet bandelet representation of an input n-dimensional signal, processes this representa- 
tion and reconstructs an output n-dimensional signal 

FIG. 2 shows, in block diagram form, an exemplary configuration of a warped wavelet packet bandelet 
processor. 

JO FIG. 3 shows, in block diagram form, an exemplary configuration of an inverse warped wavelet packet 
bandelet processor. 

FIG. 4 illustrates the partition of an image support in dyadic square regions and its quad-tree represen- 
tation. 

FIG. 5 shows, in block diagram form, an exemplary configuration of a warped wavelet packet binary 
15 tree that implements a warped wavelet packet processor. 

FIG. 6 shows, in block diagram form, a warped wavelet packet binary tree that implements a warped 
wavelet processor for n - 2 dimensional signals, over 3 levels. 

FIG, 7 shows, in block diagram form, a warped wavelet packet binary tree that implements a warped 
wavelet processor for n - 3 dimensional signals, over 2 levels. 
20 FIG. 8 shows, in block diagram form, an exemplary configuration of a Warped Subband Processor 
(WSP) along a direction d 

FIG. 9 shows, in block diagram form, an exemplary configuration of an Inverse Warped Subband Pro- 
cessor (IWSP) along direction d 

FIG. 10 shows, in block diagram form, an exemplary configuration of an inverse warped wavelet packet 
25 processor defined over the same binary tree as in FIG. 5. 

FIG. 11 shows, in block diagram form, an exemplary configuration of wavelet packet warping grid 
subsampling along the same binary tree as in FIG. 5. 

FIG. 12 shows, in block diagram form, an exemplary configuration of a bandeletisation processor. 

FIG. 13 shows, in block diagram form, an exemplary configuration of an inverse bandeletisation proces- 

30 sor, 

FIG. 14 shows, in block diagram form, an exemplary configuration of a processor unit for signal t»m- 
pression. 
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FIG. 15 shows, in block diagram form, an exemplary configuration of a processor unit for signal decom- 
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FIG. 16 shows, in block diagram form, an exemplary configuration of a processor unit for signal feature 
extraction. 

FIG. 1 7 shows, in block diagram form, an exemplary configuration of a processor unit for signal feature 
extraction computed from compressed data. 



DETAILED DESCRIPTION 

FIG. 1 shows a system exemplifying the present invention. The system takes in input an n-dimensional 
digitized signal, where n is a positive integer. The digitized signal is given by an n-dimensional array of 

10 sample values. Electrocardiograms (ECGs) are examples of a 1-dimensional signal, images are examples 
of 2-dimensional signals, and video image sequences are examples of S^limensional signals. The warped 
wavelet packet bandelet processor (101) computes from an input signal a warping geometry adapted to the 
inpufsignal, as well as bandelet coefficients that are inner products with a farniry of bandelet vectors adapted 
to the signal geometry. The processor (102) transforms the bandelet coefficients and the warping geometry 

15 according to a particular application such as signal restoration. In a transmission system, the processor 
.nay compress and code the warping geometry together with the bandelet coefficients into a bitstream, 
which is transmitted along a transmission medium, and decoded, to output processed bandelet coefficients 
and a processed geometry. The inverse warped wavelet packet bandelet processor (103) takes in input the 
processed warping geometry and the processed bandelet coefficients and computes an output signal. 

20 FIG.2 S hows,mblockdiagramfonn,^ 

processor (101). A geometrical segmentation module (201) computes from the input signal a warping 
geometry which includes parameters from which the geometrical sampling module (202) computes a list 
of signal warping grids. A signal warping grid is an array of points that are distributed according to the 
geometrical signal properties in a given region. The signal warping processor (203) computes from the input 

25 signal and the signal warping grids a list of warped coefficients that stores the signal values at the locations 
of all points of all warping grids. The output warped signalis composed of the warped coefficients together 
with the signal warping grids. The warped wavelet packet processor (204) computes a warped wavelet 
packet transform that includes warped wavelet packet coefficients and wavelet packet warping grids. The 
bandeletisation processor (205) decollates the input warped wavelet packet coefficients by applying one- 

30 dimensional linear invertible operators along selected directions of the warped wavelet packet coefficient 
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FIG. 3 shows, in block diagram form, an exemplary configuration of the inverse warped wavelet packet 
bandelet processor (103). The geometrical sampling (301) computes signal warping grids from an input 
warping geometry. The wavelet packet warping grid subsampling (302) computes wavelet packet warping 
grids from the signal warping grids. The inverse bandeletisation (303) implements the inverse of the ban- 
5 deletisation (205). It inputs the wavelet packet warping grids and bandelet coefficients and outputs a warped 
wavelet packet transform including warped wavelet packet coefficients and the corresponding wavelet packet 
waiping grids. The inverse warped wavelet packet processing applied by processor (304) implements the 
inverse of the warped wavelet packet processor (204): from the input wavelet packet transform, it computes 
a warped signal including warped coefficients and the signal warping grids. The inverse signal warping 
10 (305) recovers an output signal defined over the same sampling grid as the input signal in FIG. 2, with an 
interpolation procedure. 

In an exemplary embodiment of the present invention, each process can be carried out by a computer 
program. The computer program is run on a computing device which may comprise a storage device such 
as a magnetic disk or an optical disk on which the n-dimensional digital input signal may be stored The 

15 n-dimensional digital input signal may also have been received by the computing device after a transmission 
via a network connection such as an ethernet link, a phone line, a wireless transmission, etc... The computer 
program can be stored on the storage device or on a Read Only Memory (ROM). Each process is executed 
by a Central Processing Unit (CPU) coupled to a Random Access Memory (RAM). The output of each 
process can be stored on the storage device or transmitted via a network connection. In another exemplary 

20 embodiment, each process can be carried out by a dedicated electronic device comprising die instructions to 



Warping Grids and Warping Geometry 

The input n-dimensional discrete signal in FIG. 2 is specified by an array of samples S(m) where m = 
(mi, mn) is an n-dimensional index parameter and each m* is an integer for 1 < d < n. The position 
of these samples is represented by the point of coordinates m = (mi, mn) in W . The signal support Ss 
in 5P 1 is the halo of all m € 2" for which the signal value is defined. The halo(A'f) of a subset M of Z n is 
the subset of K" defined by a union of hypercubes; 

halo(*f)= (J [ini-l/2 l m l + l/2]x--x[m„-l/2 l m ll + l/2] 
To manipulate multidimensional indexes of arrays, the following notations are used For any integer a 
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and any positive integer b, there exists a unique pair (g, r) of a quotient and a remainder such that 

a = bq + r 

where r € {0, . . . , b — 1} (Euclidean division). We define the 'mod* and 'div' operators by a mod b = r 
and a div 6 = g. For any fc = (fci , . . . , ) € Z n , we write: 

fc modj 2 = k d mod 2 
fc divj 2 = (fei 1 ...,fcd-x>fcd div 2 l fc<M-i,...,fc f 0 

' 2<tk = (fci,..',fcd-l»2fc<f, 

fc/2<j = (fci,.. M fcd~i,fcd/2,fc<f+i,...,fcn) 

and Id = (0, • . . , 0, 1, 0, . . , , 0) is a unit displacement along direction d: the coordinate d is equal to 1 and 
all others are equal to 0. For any x € W 1 and r 6 (0, +oo) we also write 




for r < oo 



When x and u are two vectors of R™, for some m > 0, we write x.u = Y^±=i x d^d the inner product 
between two vectors x and ti. When A is a matrix of any size of entries ay, A T is the conjugate transpose 
matrix, of entries aj<, and A~ T is a shorter notation for (A T )~ l . 

The warped wavelet packet transform is computed along a union of sampling grids called warping grids. 
5 The union of the warping grids is typically not equal to the signal sampling grid and is adapted to the 
geometrical structure of the signal In an exemplary embodiment, each warping grid is represented by an 
array of points WG{%, fc) € W where i is a fixed index and fc € Z n is a position index. For a fixed region 
index t, the warping grid WG(i, fc) is only defined for a subset of values of fc, and corresponds to a region 7t< 
of the input signal support, over which the input signal has a homogeneous geometrical structure. Any other 
10 state of the art representation of arrays of sampling points may be used The warping geometry includes 
parameters calculated from the n-dimensional input signal by the geometrical segmentation module (201), 
that characterize the warping grids which are computed by the geometrical sampling module (202) from the 
warping grids. 

In an exemplary implementation, the warping geometry specifies a partition of S$ into regions Hi C BP , . 
15 with Ss = Ui and specifies a geometrical deformation 

= (<M*)> -,*>(*)) 6 ^ 
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whose support is the set of fc .€ Z n such that <Sj(ft) € Tij, For each region Hi, the geometrical sampling 
module (202) computes a signal warping grid 

Wa W )-/*<»»** l) «*. (1) 
^ nil otherwise . 

The warping geometry also specifies a regularity descriptor RegWG (t) which is a vector (pi,...,Pn) 
that characterizes the regularity of the signal samples along the different directions of the waipmg grid 
5 WG(i, fc). Along each direction d, if Pd = 0 then for fc fixed, WG(i> fc + J Id) can have irregular variations 
as a function of /. If Pd 0 then WGfa fc + Ip* I4) has regular variations as a function of I. Thus, pd = 1 
indicates that the function fc + / Id) of Z is slowly varying, whereas pd > 1 indicates that the function 

WC(i, fc + / Id) of J is nearly periodic of period pd* The vector RegWG is stored together with the warp- 
ing grids WG and is used by the bandeletisation processor (205) to perform an appropriate decorrelation 

10 of waiped wavelet packet coefficients. The geometrical segmentation module represents the partition into 
regions Hi with region parameters, die geometrical deformations £,*(fc) with deformation parameters, and 
the vector Reg W6(t) with vector parameters, from which Hi, di(k) and RegWG(i) can be recovered To 
compute the warping grids with (1), the geometrical sampling module (202) first recovers for each t, Hu 
Si(k) and RegWG (i) by inverting these operators. The following gives exemplary embodiments to com- 

15 pute the warping geometry from an n-dimensional signal These embodiments are not limitative and any 
state of the art technique may be used to compute a warping geometry from which one can compute warping 
grids. 

In an exemplary embodiment, each region Hi is the halo in IP 1 of Hi fl Z n , and is specified by a 
binary membership map for all fc € Z n bi(k) = 1 if fc € Hi and 6<(fc) = 0 otherwise. Any state of the 

20 art parameterization technique may be used to represent these binary maps, including a chain coding of 
their boundary. In yet another exemplary embodiment, a segmentation is calculated with a parameterized 
modification such as a linear warping in R n of a predefined segmentation, in order to adapt the segmentation 
to the signal The region parameters then include only the modification parameters such as the parameters 
of the linear warping operator. For images of human feces, a state of the art technique may be used to detect 

25 features such as the eyes. We then compute the affine warping which renormalizes the position of these 
features and use this affine warning to adapt a predefined segmentation of the nice into regions adapted 
to the geometrical face structures. .In yet another exemplary embodiment of the geometrical segmentation 
(201), each region Hi is a parametrized geometrical shape such as a hyperrectangle or an ellipsoid in R n 
in which case the deformation parameters can include their center (ci, c„), and their widths (w\, Wn) 

30 along each of the n direction, or Hi may be a union of such geometrical shapes. 
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In an exemplary embodiment, the signal support is divided into hyperrectangle regions K{ that are 
dyadic hypercubes of width equal to w V with j < 0, where w is a fixed integer parameter. In an exemplary 
embodiment, this partition is computed with averaged orientation vector with a splitting procedure. The 
signal gradient is a vector V5(m) = (A 2 5(mi t rri2),..., A n S(mi,Tn2)) computed at each point m € 
5 Z n nS5with: 

A d S(m) = 5(m) - S(m - l d ) for 1 < d < n. 
For any ti € R** we define the average gradient energy over a region M, in the direction of u by: 

Om(«)= £ |V5(m).«| 2 . (2) 
The average orientation is defined as the unit vector in E n that minimizes C(u): 

tiM = arg min C(u) . 
IHI=i 

We calculate the n by n matrix Q M whose coefficients are the averages for m € M of the coefficients of 
the n by n matrix (?mH = VS(m) V5(m) r . Since C(u) = v^^ti, the vector is a unit norm 
10 eigenvector of Q M corresponding to the minimum eigenvalue which is equal to Uj^ m um = Qm(«-m). 
This value measures the signal regularity in the most regular direction over M. 

The recursive splitting procedure Split(H) is defined over any bypercube 71 and initially applied to 
the whole signal support Hi =■ 5$. For any i, Split(7li) computes the average orientation ia< in M< = 
Hi n Z n . Let Ti be a parameter and d = u^Qmi^H- If C* > Ti then the square region 7^ is divided 
15 into a partition of 2 n hypercubes of twice smaller width H^u ^2»*+i-« *nd ^2»H2--ii arid SpHt(Hi) 
recursively computes SpHt(R&ii).» Split{Jlz n i^-i)* If C< < T\ then SpHt(Ki) returns the region 
index t, which is die index of the leaf corresponding to 72^ in a tree representation of the partition. In 
this tree, called 2 n -tree, any inside node has 2 n children. The region parameters of the warping geometry 
parameterizes the tree corresponding to the signal support segmentation with a state of the art representation 
20 of a 2 n -tree. In yet another exemplary embodiment, the region partition is further modified with a merging 
procedure that yield a segmentation into regions that are unions of hypercubes. Adjacent hypercube regions 
are merged using a criterion similar to the splitting criterion above: whenever the minimum eigenvalue of 
Q M is smaller than a fixed threshold T 2 , where At is the union of two region grids M{ x = %i 1 n Z n and 
Mi 2 = %i 7 H Z n , the two adjacent regions i\ and i 2 are merged The region parameters then includes 
25 parameters specifying this merging. In dimension n = 2, the partition is represented by a quad-tree, as 
illustrated in FIG. 4. The square image support is subdivided into four equal square regions. The lower 
left square is not subdivided and corresponds to the leaf 4 of the quad-tree above. The upper left square is 
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subdivided in four squares corresponding to the leaves 20, 21, 22, 23 of the quad-tree. In an example of a 
region merging procedure, regions 20 and 21 are meiged together into a single region, and regions 6, 30 and 
126 are also merged together 

In dimension n = 3, in yet another exemplary embodiment of the geometrical segmentation (201), the 
5 regions %i are 3-dimensional rectangles, of fixed width w along one of the direction, say z 3 , and whose 
projections on the plane (ari,X2) are dyadic squares that define a partition of the signal support which 
is represented by a quad-tree. This quad- tree can be calculated with a splitting procedure similar, to the 
quad-tree splitting procedure for n = 2, using average gradient energy measurements over each of the 
3-dimensional rectangles specified by a square in the plane (xi , X2) and a width w along X3. 

10 Directional geometrical deformations are exemplary embodiments of geometrical deformations 8 = 
Ww (*)* **»n (*)) tnat can flave * + 1 different directions. The direction 0 corresponds to a grid parallel to 
Z n for which 5,- (fc) = k + a, where c,- is a constant vector. For any 1 < d < n, a geometrical deformation 
along direction d satisfies S it i(k) = k ( -h for / ^ d and S i/t {k) = k* + c^fci, k d -i , kj+u fen)- In 
an exemplary embodiment, the constants cy for J ^ d are chosen to depend only upon the geometry of 

15 Tsadrregion**?^ and their value do not need to be stored It may be minus the minimum coordinate in the 
direction I among all points in Hi. Any state of the art parameterization technique may be used to define the 
deformation parameters that specify the function aj. In an exemplary embodiment, deformation parameters 
are decomposition coefficients of c^ffa,. &<m-i> in a basis such as a wavelet basis or a block 
cosine basis defined over the set of parameter indexes (k\ , kd-i , Atf+i , *n) fr° m which there exists fcd 

20 such that (ki + c* t i, kd-i + Ctf-u Ad , A^+i + c^i, kn + q |fl ) € 7£<. Fast transforms as described 
in S. Mallat, "A wavelet tour of signal processing", Academic Press, 1998, can be used to compute these 
decomposition coefficients. 

In an exemplary embodiment, for each regions Hi of a signal partition, the deformation direction is 
computed by finding the direction in which the signal is the most irregular. It is done by calculating the 

25 matrix ^ previously defined. If the trace of is smaller man an adjustable parameter T* then the 
region is uniformly regular, in which case the geometrical deformation has the direction 0 and RegWG(t) = 
(1, 1, 1). Otherwise the direction in which the signal is the most irregular is defined as the direction d 
such mat 

The regularity descriptor RegWG(t) = (pi,P2, —>Pn) is then defined by pi « 1 for each / ^ d and pa = 0. 
30 We shall see that the value pa — 0 can later be updated if there exist quasi-periodic variations along the 
direction <L 
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In dimension n = 2, a directional geometrical deformation in the direction 1 can be written 

**(*)» (*1 + (fc)»*2 + Ctf) ' 

and in the direction 2 

**(*) = (to +<*i ,fc2 + C|^(fci))- 

In an exemplary embodiment, the geometrical deformation over each region is also calculated from average 
orientation vectors. Let us consider a region whose deformation direction is d = 2, the other case d = 1 
5 being treated similarly. For each J € Z let be the set of integers 7712 such that (I, m2) € ft*. The first 
parameter of the geometrical deformation c^i is an integer that is set to 0. To compute C{j(ki), we compute 
the average orientation vector U{(1) = (1*1,1 (0»H2(0) ovcr ^M- Let o be the rninimum of ail integers I 
such mat fcy # 0. We define 

The discrete signal 6^2 (mi) is then decomposed over a basis such as a discrete wavelet basis or a discrete 
10 cosine basis with a linear operator. To obtain a sparse parameterization, in an exemplary embodiment, all 
decomposition coefficients of c^mi) of amplitude below a threshold T 8 are set to zero. For compression 
applications, in yet another exemplary embodiment, the decomposition coefficients of c^mi) are quan- 
tized by a scalar quantization. We define Ct^(mi) as the signal reconstructed by applying the inverse linear 
operator on the thresholded or quantized coefficients of ^(mx). The deformation parameters of the warp- 
15 ing geometry specify the direction d m 2 and the thresholded or quantized decomposition coefficients used 
to compute (^(fci). An equivalent procedure computes the geometrical deformation in the direction d = 1. 

In yet another exemplary embodiment, in any region Tlx C K", a directional geometrical deformation 
is computed in dimension n in a direction 1 < d • < n by detecting local maxima of the gradient The 
signal gradient V5(Jfe) is calculated at each h e Z n R £5. For any p ^ d we set S^k) = Jbp. For 
20 any fixed (fci,...,fc*-i,*a+i,^ € Tti 

and if this set is not empty, we find in this set the integer &<* such that | V*S(fc)| 2 is maximum in k = 
(tot kd-u kdi > *n) 211(3 define 

We decompose c,*^ in a discrete basis such as a wavelet basis in dimension n - 1 and threshold its coefficients 
or quantize its coefficients and denote by c^Afi, kd-i , fcd+i , fcr») the hypersurface reconstructed from 
25 the thresholded or quantized coefficients of c,- f< f . The deformation parameters specify the direction d and the 
thresholded or quantized decomposition coefficients used to compute c^ki). 
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For an input video image sequence, which is an n = 3 dimensional signal, in yet another exemplary 
embodiment, the warping grid segmentation module (201) computes a directional geometric deformation by 
applying a time displacement to a 2-dimensional geometrical defonnation using estimated motion vectors. 
The video sequence is written S(mi } m2,mz) where m$ is die time parameter, and for a fixed n%z = J, 
5 S(mx, m 2) /) is an image. For any region TU and m 3 = / fixed, let TJy be the image region of all (mi,m 2 ) 
such that (mi, ma, /) € H%. To a three-dimensional region 7^, we associate a spatial geometrical defor- 
mation (ci f i(*i . k 2 ) 5 c;o(k u fco))> and a two-dimensional regularity descriptor (pi,p2) which indicates the 
spatial direction in which the signal is regular. Let a be the minimum I for which 71^ ^0. In an exemplary 
embodiment, the said spatial geometrical deformation and two dimensional regularity descriptor (pi,P2) are 

10 provided by computing a two-dimensional directional geometrical deformation in die slice 71^. We define 
RegWG = (pi,J>2» 1) which indicates that the variations in time are regular. For any J for which Kij is 
not empty, an average motion vector is computed over TL^ by calculating the average displacement of the 
signal components in this region with respect to their position in die previous image S(mi, m*,/ - 1). Any 
state of the art motion estimation algorithm may be used, including matching procedures or gradient based 

15 computations. Let Vi{l) — (v<,i(/),t>t,2(0) be the average motion vector of the region and 

In an exemplary embodiment, a sparse parameterization is obtained by decomposing this displacement vec- 
tor in a basis, thresholding or quantizing the decomposition coefficients and reconstructing an approximate 
displacement vector (^,1(1713), 1^(1713)) from the quantized or thresholded coefficients. Defonnation 
parameters include these thresholded or quantized decomposition coefficients. The motion geometric defor- 
20 mation is then defined by 

Si{k) = (c,i(*i, *2) + Wi A (k z ) ,c^(fci, k 2 ) + Wij(h) , fc 3 ) • 

The geometrical segmentation module (201) can also detect the existence of quasi-periodic signal varia- 
tions, measures these quasi-periods to adapt the warping grid in each region. We shall first consider tie case 
of an input one-dimensional signal and then describe an exemplary embodiment for n-dimensional signals. 
For a one-dimensional signal 5(m), the geometrical segmentation processor divides the signal support 
25 Ss into intervals [<h , a*+i] in which either the signal is quasi-periodic, with quasi-periods that have a 
bounded relative variation, or in which the signal has no quasi-periodicity. The value of RegWG (t) is used 
to indicate what is the nature of the signal in each interval. RegWG (i) = 0 indicates that S is not quasi- 
periodic over [ai , a< + i - 1], in which case the geometric deformation mapping is simply <*{(k) = * + a for 
some integer c,-. RegWG (i) = 2 1 for some integer / indicates that 5 is quasi-periodic inside [a* , Of +1 - 1] 
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the varying period being close to 2*. In a particular embodiment, the time-varying period in [a* , a,+i - 1] 
is defined from J\ consecutive period measurements Pi(j) for 1 < j < J { with J2jLi p iU) - ^+1 ~ 
1 - Of, To obtain a warping grid WG(i 1 k) such that the corresponding waiped signal WC(i, k) has a 
quasi-period equal to 2', the geometrical deformation 6{(k) is chosen to be a mapping from an interval 
5 [0, Ji2 f ] to [ai , Ot+i - 1] that maps 0 to a* and each j2 l to a, + P*(l) + • • • + P*(j). Let 6q = <H and 
bj = a, + P»(1) + *•■ + -P.t?) for all j in {1,.. M J<}. We set 8i{j$) = for; = 0,...,Ji and 
intermediate values 5i{k) fork ^ are calculated with any state of the art interpolation procedure such as 
the piecewise linear one: 

Vk€b*2 , + l,0 , + l)2 l -l] , 6 i (k)^b j + {b j ^ 1 -b j )x(2' l k-j) . 

Any state of the art periodic detection procedure may be used to segment the signal support and compute 
10 quasi-periods Pi(j) ^ interval. In an exemplary embodiment, the periodicity is measured with a 
variogram, defined over [0, L] and for any (m,p) € Z 2 : 

^Kp) = £ M m + 0 - 5(m + p + 1)| . (3) 

; - ~ • i=q - — ..... 

Period values are estimated in a predefined range [Pmfa, Pym)* An admissible period computed from ref- 
erence point m is defined as the minimum p > P m ^ with p < Pmax such that V(m,p) < T where T is a 
fixed threshold. If such a j? does not exist then S(m) is considered as non-periodic on [m, m + Pmax]. 
i5 We suppose that the signal support Ss is an interval and cut iteratively this interval into regions 
[a,-, Oj + i - 1], by constructing iteratively the sequence (a*)*. The sequence is initialized by setting oi to 
be die first point of Ss* At each iteration, a value a< has been determined and we compute the next value 
c4+i as follows. 

• If the variogram (3) indicates that S is not periodic on [a*, a* + Pmax]» i-e. there is no admissible 
20 period from reference point m = a,-, then RegWG(t) = 0, We find the minimum integer q for which 

when m = a* + qPmax + 1 the variogram (3) finds an admissible period p from reference point 
rn = Oi + q Pmax + 1 and we set a* + i = m - 1. 

• If the variogram (3) finds an admissible period p computed from reference point m = o< then we set 
RegWGft) a* 2 l where / is die closest integer to log 2 (p). We set jft(l) = p, to = <H» &i = <H + P*(l) 

25 and compute iteratively die next values of b q , for q > 1 as follows. We suppose that b q has been 

computed. Every time the variogram (3) finds an admissible period p from reference point m = b q> 
and if max(2 l ,p)/xnin(2 i ,p) < T 4 where T 4 is a fixed threshold then we set P<(g + 1) = p and 
ss 6 ? + p. This is continued until the variogram (3) does not find an admissible period p from 
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some reference point m = bj { or if it finds a period p such that max(2 l ,p) / min(2 i ,p) > T 4 . The 
quasi-period sequence P{(j) for 1 < k < K{ is considered as a signal that decomposed in a basis 
and whose decomposition coefficients are thresholded or quantized We write Pi(j) the sequence 
reconstructed from these thresholded or quantized decomposition coefficients. We then set a^i = 
5 Oi + ZlUPiV)- 

This progressive construction of intervals [ai, Of+i — 1] continues until we reach the end of the support $$. 
In an exemplary embodiment, the region and deformation parameters of the warping geometry specify for 
each interval the size 2 1 , the number of quasi-periods J<, and the thresholded or quantized decomposition 
coefficients used to compute the quasi-period sequence Pi{j)- When the signal is not periodic in [a^, Of+i — ], 
10 the length Pf (1) - a* + i - a* is stored. Since a^+i = a* + J2jL\ flfO") interval boundaries are recovered 
from these parameters. 

We now consider for n > 1 an n-dimensional signal S(m) whose support Ss has been partitioned 
into regions Hi over which a geometrical deformation £<(fc) was calculated together with RegWG(t) = 
(pii —iPn) with pi = 1 or pi = 0. Several exemplary embodiments were described to compute such a 
15 partition. When the correlation type RegWG(i) of a region i is such thatp* = 6 for a single direction d, an 
exemplary embodiment of the geometrical segmentation module (201) scans the signal in region i along the 
direction d for which pd = 0, and if there exists some quasi-periodicity along the samples of the warping 
grid defined by ft farther divides %i according to these quasi-periodicity. The quasi-periodicity 

is measured with a variogram in the direction d t with a similar strategy as the exemplary embodiment 
20 described for n = 1 dimensional signals. This leads to new sub-regions inside each of which a new 
geometric deformation further modifies the geometrical deformation 5,(fc) so that the resulting warping 
grid has a quasi-periodicity of fixed size 2 1 along the direction d, and we replace die <f-th component of 
RegWG(t, j) = (pi, ...,Pn) with pd = 2 l . When the number of directions d such that pa = 0 is higher 
than 1, a similar periodicity scanning and segmentation is performed, wherein the periodicity scanning 
25 from a reference point consists in finding a family of independent period vectors by minimizing a variogram 
function, and wherein the periodicity scanning is performed along the samples of the warping grid defined by 
6i(k) along the axes for which pd = 0. After proper modification of the geometrical deformation function 
tf(i,fc), the signal is quasi-periodic of period 2^ along some directions d of the warping grid defined by 
5 (i, fc), and for these directions the vector RegWG(s) is updated by replacing its d-th coordinate previously 
30 equal toO with 2 / *. 
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Signal warping and inverse warping 

The signal warping processor (203) computes a value WC(i, k) of the signal at each warping grid points 
WG(i, k) ^ nil, given the input signal sample values S(m). This calculation can be performed with any 
state of the art interpolation procedure. The output warped signal includes the warped coefficients WC(i t k) 
5 together with the signal warping grid WG{%> fc). 

An exemplary embodiment is obtained by computing the interpolation with an interpolation function 
4>{x) for x e IP 1 which satisfies #0) = 1 and <f>(m) = 0 for m € Z n and m ^ 0. At any location 
WG{%) k) — x^ nil, the interpolated coefficients are then 

WC(i t k) = £ S(m)fl«-m). (4) 

An example of an interpolation function <f>{x) is a separable function 

<f>(x u ...,x„) = ^o(a?i) ... <fo(z n ) , (5) 

10 where ^o(<) is an interpolation function with <f>o(0) = 1 and ^o(m) = 0 for m e Z and m ^ 0. In an 
exemplary embodiment, an interpolation procedure such as described in P. Thevenaz, T. Blu, M. Unser, 
"Interpolation Revisited," IEEE Transactions on Medical Imaging, vol 19, no. 7, pp. 739-758, July 2000, 
is used to compute WC(i, &). State of the art techniques are used to adapt the interpolation kernel at the 
boundary of Ss< 

15 The inverse signal warping processor (305) in FIG. 3 takes in input a warped signal including the 
signal warping grids WG\ and warped coefficients WC\ at the corresponding locations, and computes fee 
output signal that is written S 0 {m) defined for all m € Z n in the support Ss of the input signal in FIG. 
2. The operation it performs is an approximation of the inverse of the interpolation implemented by the 
signal warping processor (203). The output signal can be calculated with any state of the art interpolation 
20 procedure. An exemplary embodiment of this interpolation is calculated as follows. For each m € Z n C\Ss, 
we find the warping grid point WG(i mi k m ) = x that has a minimum distance d(x, m) where d may be 
any distance such as the Euclidean distance. The value S c (m) is computed by finding local coordinates c 
of m in the warping grid WG(i, k) and performing an interpolation with a uniform metric in this warping 
grid. In an exemplary embodiment, for each direction 1 < d < n, we find among the two neighbors 
25 WG{im,k m ± l d ) the one WG(im, k m + c^U) with a = 1 or -1 that is the closest to m. The local 
coordinates c = fa, Cn) are then defined as the sums of tire coordinates of fc m and the coordinates of 
the vector m - WG(i m , km) in the basis of the n vectors e d (W r G(i m , + c d l d ) - WG^ k m )) for 
> 1 < d < n. 
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Given the local coordinate c of the point m 6 Z n in the warping grid WG(i m , k) 9 in an exemplary 
embodiment, S 0 (m) is computed with an interpolation function <f>(z) which satisfies ^(0) = 1 and ^(fc) = 
0 if k € Z n . The interpolated coefficients are then given by 

In a preferred embodiment, the interpolation function <f>(x) is a separable function like in (5). State of the 
5 art techniques are used to adapt the interpolation kernel at the boundary of the waiping grid WG(i m , fe). 

Warped Wavelet Packet Processor and Inverse 

The warped wavelet packet processor (204) computes a warped wavelet packet transform from the input 
warped signal, with a binary tree of one-dimensional warped subband processings. An exemplary embod- 
iment of this binary tree is illustrated in FIG. 5. The input warped signal includes the warped coefficients 

10 WCx and the signal warping grids WG\. Similarly to a tree-structured subband coder as in MJ. Smith and 
T.P. Barnwell ID, **Exact reconstruction for tree-structured subband coders", IEEE Trans. Acoust, Speech 
and Signal Proa, voL 34, no. 3, pp. 431-441, 1986, a warped wavelet packet transform is calculated with 
a binary tree of subband decompositions which splits the input signal into two signals which have the same 
total number of samples as the input signal, but in this case the subband decomposition is performed by a 

15 Warped Subband Processor (WSP) adapted to the warping grids that are typically different from the input 
signal sampling grid This warped subband processor thus operates differently from a state of the art sub- 
band filtering procedure, as in S. Mallat, "A theory for multiresolution signal decomposition : the wavelet 
representation," IEEE Transaction on Pattern Analysis and Machine Intelligence, vol. 1 1, p. 674-693, July 
1989. Each node of the binary tree is labeled by an integer p. At the root p = 1. The left and right children 

20 of a node p are respectively labeled by the integers 2p and 2p+ 1. Each Warped Subband Processor (WSP) in 
(501), (502), (503), (504) and (505) performs a one-dimensional warped subband processing along the lines 
of the input warped coefficients in a direction 1 < dp < n, which can be chosen arbitrarily for each node p 
of the binary tree. At a node p> this warped subband processing splits the input warping grids WG P and the 
corresponding warped coefficients WC P into two sets of warping grids WGzp and W(?2p+i and two sets 

25 of warped coefficients WC^p and WQtp+\ defined over these grids. The warped coefficients WC^p are ob- 
tained with a warped subband filtering using a low-pass filter whereas the warped coefficients WC2p+i are 
obtained with a high-pass filter. The output warped wavelet packet transform corresponds to the family of 
warped coefficients WC P together with their warping grids WG pt at the nodes p € P which are the leaves 
of the warped wavelet packet binary tree, In FIG. 5, it is composed of WC* and WG^ WCj and WG 7 , 
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WCio and WGio, WCu and WG n , WC X2 and WG U , WC& and WG XZ > so P = {4, 7, 10, 11, 12, 13} 
in this case. 

A warped wavelet processor is a particular embodiment of waiped wavelet packet processor, corre- 
sponding to a specific binary tree. This binary tree is the same as the one used to compute a separable 
5 n-dimensional non-warped wavelet transform, as in S. Mallat, "A wavelet tour of signal processing", Aca- 
demic Press, 1998. In dimension n = 2, the binary tree is illustrated in FIG. 6. This binary tree is composed 
of WSP (Warped Subband Processor) (601-609) in directions d p = 1 if 2* < p < 2* +1 with k even, and 
d p z=2 if 2* < p < 2* +1 with k odd After the first two levels of WSP decompositions in the binary tree, 
only the warped coefficients WC4 and waiping grids WG4 are further processed The same succession of 

10 warped subband processors in the directions d = 1 and d = 2 is then applied and only the output of the 
low-pass branch corresponding to WC P and WG P for p = 4* is further sub-decomposed, as illustrated in 
FIG. 6. This binary tree of processors can be further continued on more levels. 

In dimension n = 3, the binary tree of a warped wavelet processor is illustrated in FIG. 7. This binary 
tree is composed of WSP (701-714) in directions dp = 1 if 2* < p < 2 fc+1 with k « 0 mod 3, dp = 2 

15 if 2* <p< 2 fc+1 with k «= 1 mod 3, dp = 3 if 2* < p < 2* +1 with fc = 2 mod 3. After three levels 
of WSP decompositions only the warped coefficients WC P and warping grids WG P with p = 8 fc , obtained 
wifli a succession of low-pass subband waiped filtering, are further processed, like WC$ and WG* in FIG. 
7. 

An exemplary configuration of a warped subband processor in a direction d is illustrated by the block 
20 diagram in FIG. 8. Hie grid subsampling along direction d performed by module (801) splits the warping 
grid WG P into two warping grids WGzp and W(?2jh-i» defining for all % and all fc 

WG2p(hk) = WG p (i,2 d k) and WG 2p +i(i,k) = WG p {i,2 d k + l d ) , 

and ' 
RegWG2p(t) = IUgWG w (i) - RegWG p (i)A, . 

Note that in this last formula, the resulting vectors RegWG ^(i) and RegWG^xW may have non integer 
entries, like 1/2 or 1/4. In the description of the present invention, waiping grids WG P corresponding to 
separate nodes of the wavelet packet binary tree are considered as stored in separate arrays, to simplify the 
explanations, but can also be stored in a single array. 

If the warped wavelet packet transform is done individually within each grid, when reconstructing an 
output signal from processed warped coefficients, it creates discontinuities at the grid boundaries. To avoid 
this effect, the warped wavelet packet processor performs the subband filtering across warping grids. Grid 
connections are computed by the neighborhood connection processor (802), which establishes neighborhood 
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connections between boundary points of different regions along the direction <L A warped subband filtering 
processor (803) along direction d takes in input warped wavelet coefficients WC P with the corresponding 
warping grids WG P and a neighborhood connection graph WN p (i, fc, A, d) along direction d and performs 
a subband filtering within each grid of coefficients WC p (i y k) and across grids, with adapted filters. It 
5 outputs the subband warped coefficients WC^ and WC^\. 

Hie neighborhood connection processor (802) establishes a neighborhood relation which is symmetric, 
meaning that if a first point is a left neighbor of a second point then the second point is the right neighbor 
of the first point In an exemplary embodiment, for each WG p (i,k) ^ nil, WN p (i,k 9 -l,d) = (ij,fcj) 
gives the index *"/ of the region where the left neighbor of WG p (i, fc) in the direction d is located and fcj 

10 is its index in this region. Similarly, WN p (i>k, l,d) = (i r ,fc r ) gives the index i r of the region where the 
right neighbor of (i, fc) in the direction d is located and k r is its index in this region. The computation 
of WN p {i,kt±l,d) proceeds as follows: if WG p {i,k - l d ) ^ nil then WN^k.-l.d) is set equal to 
(i, fc - l<f ), and otherwise (i, fc) is included in the Right list that stores the indexes of all points that do not 
have a left neighbor in the same region. Similarly, if WG p {i y k + l d ) # nil then W\Np(t\fc,l,d) is set 

15 equal to («, fc + lrf), and otherwise (i, fc) is included in the Left list that includes the indexes of all points 
that do not have a right neighbor in the same region. 

For all (i, fc) in the Right list we search for a left neighbor point (*', fc') in the Left list, with i' ^ % 
and which minimizes Coat^VG p (i\ fc'), WG p (i, fc)), where Cost (a, z') is a cost function defined over 
paiis of points in RP. We then set WN p (i,k,-l,d) = (i',k'). Similarly, for all the points (t,k) in 

20 the Left list we search for a right neighbor point (i'tfc') in the Right list, with t' ^ *, which minimizes 
CortiWG^hk^WGptftk')) and set WN p (i,k,l,d) = (t^tf). To guarantee that the neighborhood 
connection graph is symmetrical, Le. that 

WN p (i, k.l.d) = (f,k') * WN p (i',k\ -1,4) = (i.fc) 

for each point (i, fc) of tbe Left list, if WN p (i, fc, 1, d) = (»', fc') and WN p (i', fc', -1, d) ^ (», fc) we set 
WTVpfi, fc, 1, d) = nO. which means that it has no right neighbor. Similarly, for each point («, fc) of die 
25 Right list, if WN p (i, fc, -1, d) = (?, fc') and ^JVp(i', fc', l,d) # (t, fc) we set WN p (i, fc, -1, d) = niL 
A first example of cost function is derived from a norm in R" for some r > 0 

Cost(WG p (i,k), WG p (i',k')) m ||FT<?,(i,fc) - WG p (i',k')\\ r . 

A second example computes a distance after translating the point in the Left list by a prediction vector 
in the direction d 

CoatiWG^klWG^.k')) = \\WG p (i,k) + n* - WG p (i\k')\\ r 
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where we may choose = WG p (i,k) - WG p (i,k - l<j). Any other cost function may be used to 
compute the neighborhood connection graph WN p . 

The warped subband filtering processor (803) performs a one-dimensional subband filtering of the 
waiped coefficients WC P along the direction d, across the different grids, through the connections WN P of 
5 the waiping grids WG P and outputs WCi p and WChp+i. The inverse aligned subband filtering processor 
(903) in FIG. 9 implements fee inverse by reconstructing WC P from WChp and WC^i, given WG P and 
WN P . The detailed implementation of both processors is feus described together. 

A state of fee art one-dimensional uniform subband filtering is computed with a Finite Impulse Response 
low-pass filter B and an FIR high pass filter (?, which admit a dual pair of FIR reconstruction filters H and 
10 (?, as described in S. Mallat, "A wavelet tour of signal processing", Academic Press, 1998. The subband 
decomposition performs a filtering with H and G of WC p {i, k) along fee direction d and subsamples fee 
output on even index points within each grid. 

Let Sh and Sg be the finite support of the filters E(l) and G(l) for / 6 Z. If (i, fe) satisfies 

VI e S H and Wl 6 S G WG p {i,2 d k + ll d ) ? nil (6) 
then a state of the art uniform subband filtering along fee direction d computes: 

15 and 

WChp+iih k) = 2 <?(') fFCM*. 2d* + »d) • W 

For fe = *<f, —i *n) we write ^ + *<f m & s d + suPP° rt of JJ and G? translated by The 

reconstruction is calculated wife: 

wCp{i,k)= nw-kjwcjpfrk-iu) (9) 
+ £ 6(a-hi)wc*Hlh*-m). (io) 

For signals of dimension n > 2, a phase aligned subband filtering computes WC^p and WChp+i from 
WC P so that for all («, fc) which satisfy (6) the equality (7) and (8) are satisfied. The subband filtering 
processor (803) is then said to have a phase alignment coherent wife fee signal warping grids* To maintain 
this phase alignment property within each grid t, fee subband filtering processor (803) adapts the filtering 
20 at the interface between warping grids. The inverse warped subband filtering processor (903) computes 
fee inverse subband filtering (9) within each grid i and adapts fee inverse filtering at fee interface between 
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warping grids. For signals of dimension n = 1, the warped subband filtering processor (803) does not 
necessarily maintain the phase alignment property, which means mat local shift may be introduced when 
calculating the warped coefficients. 

In an exemplary implementation, the one-dimensional subband filtering along d is performed by the 
5 processors (803) and (903), by decomposing WG P as a union of connected lines along the direction d, as 
follows. In both processors, we traverse WN P to build an array Begin(q) = (», k) corresponding to grid 
points that do not have a left neighbor, which means that WN p (i t k, -1, d) = nil Each line indexed by g is 
initialized with L q (0) = Begin(q) if k mod,, 2 = 0 and 1,(1) = Begin{q) if k mod* 2 = 1. Then for any 
m > 0 we iteratively compute the right neighbor of X,(m) = (f , *) and set L q (m + 1) = WN p {i, k, 1, d) 
1 0 and stop when it is nil For each q fixed, we compute (7) and (8), with a direct subband filtering of L,. 

We maintain the phase alignment property by imposing that even samples of L q correspond to even 
samples of the grids WG p (i,k) along the direction d: if (t.Jfe) = L p (m) then Jfemod<,2 = mmod2. 
To enforce this property, points are inserted with the following procedure. Because of our initialization 
procedure, the first point X,(0) or i,(l) satisfies this alignment property. Then we traverse L q and insert 
15 a-new point- with-a- special-symbol Insert when the phase alignment is not respected For all m > 1, if 
(i, k) = L q (m) and k mod,, 2 ^ m mod 2 then for all / > m the entry L q (l) is shifted to £,(! + 1) and 
we set L q (m) = Insert, until the end of the line where L q (m) = ml In dimension n = 1, not maintaining 
the phase alignment property for me processors (803) and (903) means not inserting new points along each 
line L q (m), in which case we may have (i, k) = L p (m) and k modV 2^m mod 2. 
20 To implement the processor (803) the warped coefficients values WC p (i t k) along the line L q (m) are 
stored in C q (m). If L q {m) # nil and L q {m) ± Insert then C q (m) = WC p {L q {m)). If L q (m) = (», k) is a 
point sufficiently far from the grid boundary and satisfies (6) then the values C q (m) are filtered with the same 
procedure as in (7) and (8). Otherwise, the filters are adapted to the interface between the connected grids 
to guarantee the perfect reconstruction property. In an exemplary implementation, the connecting filters are 
25 implemented by inserting new samples at grid connections where L q (m) = Insert and using these inserted 
values in a filtering scheme using the same filters H and G. Each inserted value C q (m) is expressed in terms 
of all other values C q {1), including other inserted coefficients, with linear interpolation coefficients I m (l). 
Let N(q) be the number of inserted values where L q {m) = Insert along a line L q . We call S q the set of 
indexes m such mat L q (m) ? nil and S iiq the set of indexes m such mat L q {m) = Insert. The unknown 
30 value of all these inserted coefficients are calculated by solving the N(q) by N(q) linear system: 

Vm € S iA > = J2 7 -W C « W ° 0 (U) 

with I m {m) = -1. If the linear system (1 1) does not have a unique solution, then we compute the vector of 
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inserted coefficients (C q {m)) m ^ of minimum Euclidean norm which minimizes Eme^ K™)| 2 > with 
a singular value decomposition. In a preferred embodiment, J2ies q ^m(0 = 0 and the coefficients I m (l) are 
convolution kernel coefficients I m (l) = I(m - /), where border problems when defining I m for extreme 
values of m are solved with a state of the art technique. If !(-/) = J(/), one may use a a symmetric 
5 extension of C q {l) beyond its support In an exemplary embodiment we choose 1(1) = 1/2 if I = 1 and 
I s — l f and /(/) = 0 if |(| > 1. Since there is no consecutive inserted coefficients in L q (m) in this case the 
solution of (1 1) is computed explicitly by 

Warped subband coefficients are then calculated with a subband filtering of the line values: 

-Df (2m) = '(0 c *( 2m + 0 ( 12 ) 

and 

D q {2m + 1) = Y, °(0 C *( 2m + 0 • ( 13 ) 
fefio 

J0 Then values of PTC^ and WC^p+i arrays are computed as follows: if L q (2m) = (i, fc) ^ Insert then 
TFCfc,(i, * div d 2) = JD ff (2m) and if £ fl (2m+l) = (t, fc) ^ Insert then TV P Ci p+1 (t 1 fe div d 2) == D q (2m+ 
1). Observe that the coefficients D q (m) corresponding to inserted coefficients, i.e. such that L q (m) = 
Insert, are discarded If (t , &) satisfies (6) then (12) and (13) give the same result as (7) and (8). State 
of the art boundary techniques explained for example in S. Mallat, **A wavelet tour of signal processing;", 

IS Academic Press, 1 998, such as symmetric procedures are used to solve boundary issues in the convolutions 
(12) and (13). In an exemplary implementation, these convolutions are computed with a state of the art 
lifting scheme as in L Daubechies and W. Sweldens, "Factoring wavelet transforms into lifting steps", Jour, 
of Fourier Anal. AppL, Vol. 4, No. 3, pp. 245-267, 1998. In an exemplary implementation, the filters H 
and G are chosen to be the 7/9 Cohen, Daubechies, Feauveau filters specified in the above reference. 

20 For a fixed q 9 let M(q) be the size along m of the vectors C q (m) and D q (m) and N(q) be the number 
of inserted values corresponding to m 6 Sf rf . We denote by O q the linear operator from R M M to tt*^ 
which associates to any vector C q (m) the vector D q (m) with the filtering and subsampling formulas (12) 
and (13). The values of the line q in the direction d of the input WC P and of the outputs WC 2p and WCip+i 
are stored in the vectors C q {m) and Z? g (m) of size M(q) - N{q), that equal respectively to C q {m) and 

25 D q (m) for all m € S q - S^. Inserting coefficients according to (1 1), computing the coefficients D q (m) 
according to formulas (12) and (13) and then discarding coefficients D q (m) corresponding to inserted coef- 
ficients for indexes m e S iA defines an operator O q ftom <<rM%) to R"(«W%> that maps the vector 
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(C ff (m)) m€ ^-5j^ t0 CDfffaJW^-Sfrf- This linear operator O q performs an adaptive filtering, with mod- 
ified filter coefficients over the points L q (m) = (t,fc) which do not satisfy (6). Any other state of the 
art implementation of O q , including a direct lifting scheme with no coefficient insertion, may be used to 
compute D q (m) and hence each line q of WC^ and WC^jy+v 
5 The inverse processor (903), recovers the waiped coefficients values WC P along each line L q (m) from 
the values of WC^p and WC^+i along the corresponding lines. The vectors I m are defined by J m = 
O q T I m . Their expression is then for indexes m away from the boundaries: 

/m(2r)= X^(0im(2r-0 (14) 
test 

and 

4(2r + 1) = Y, *(0 Jm(2r - 1) . (15) 

When the coefficients J m (/) are convolution kernel coefficients, the coefficients of I m (l) have a simple 
10 structure: I m {l) = Iq{1 - m) for even m and I TO (/) s - m) far odd m, where J 0 and Ji are two 
filters. This is not true near-the- boundaries,- where the computation of J m = O q T I m is obtained by using 
appropriate boundary techniques when evaluating the sum (14) and (15), which are apparent to those skilled 
in die art 

Each array D q is defined for indexes m 6 S q - S iiq as follows: if L q (2m) = (t,fc) ^ Insert 
15 then D g (2m) = WChp{i,kdiv d 2) and if JS,(2m + 1) = (i,fe) 96 Insert then X> 9 (2m + 1) m 
W(hp+\{h * divj 2). To compute Dg(ro) for all m 6 5^ we solve the system of JV(g) linear equations 
with N(q) unknowns: 

, e(m) = £l ro (/)I>,(/)-0. (16) 

If die linear system (16) does not have a unique solution, then we compute the vector of inserted coef- 
ficients (Dq(m)) m e$ iiq of minimum Euclidean norm which minimizes J2mGSi q |e(ni)| 2 , with a singular 
20 value decomposition. 

An inverse subband filtering is then applied on the line values D q (m). For all / £ S q - Si A 

Cfl(0= £ B(l-2m)D q (2m)+ £ <?(/- 2m)D q (2m + 1) (17) 

and WCy(£c (m)) = C q (m). State of the art boundary techniques, as described in S. Mallat, "A wavelet 
tour of signal processing", Academic Press, 1998, such as symmetric procedures are used to solve boundary 
issues in die two convolutions (17). In an exemplary implementation, these convolutions are computed wife 
25 a state of the art lifting scheme as in I. Daubechies and W. Sweldens, "Factoring wavelet transforms into 
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lifting steps", Jour, of Fourier Anal. Appl., Vol 4, No. 3, pp. 245-267, 1998. For each q, the insertion 
system (16) together with (17) computes the inverse O^ 1 of the linear operator O q previously defined Any 
other implementation of 0~ l may be used, including a lifting scheme. In dimension n = 1 if the phase 
alignment property is not maintained then there is no inserted coefficients and the system (16) is therefore 
5 not used. 

In yet another implementation of the warped subband filtering processor (803) in direction d and its 
inverse (903), the subband filtering takes into account the fact that points in the sampling grids WG p (i,k) 
are irregularly spaced, using subband filtering on irregular point sets. In this embodiment, the position 
(abscissa) of each point of a line L,(m) is stored in * 9 (m). For all {q,m) with L q {m) ^ Insert and 
10 L q {m) ? nil, we set X q (m) = WG p {L q (m)).l d , or any curvilinear abscissa defined along the polygonal 
line joining the points (WG p (L q (m))) m es,-Si4- If «* 6 S iA then 

*»M = \ (*«("* " !) + *M + !)) 

and if X,(m) = nil then X q {m) = nfl. 

To implement the warped subband filtering processor (803), the value of C q (m) at theJV (g) inserted 
positions where £,(m) = Insert are calculated by solving an N(q) by N{q) linear system: 

Vme5^ , £/ m (i)C7,(/) = 0 (18) 

IS with I m (m) = -1. If this system does not have a single solution, the inserted values are calculated with a 
singular value decomposition as previously explained For / # m, the values of 7 m (i) may depend upon m 
and be obtained by any state of fee art interpolation kernel for irregular spaced samples located at positions 
X q (m), such as a Lagrange interpolation. 

The low-pass and high-pass convolutions (12) and (13) are then replaced by a state of fee art subband 

20 filtering procedure for irregular point sets. In an exemplary implementation, such a subband filtering, also 
called subdivision scheme, or one step wavelet transform, can be calculated with a lifting scheme, with 
any state of fee art procedure such as the ID subdivision in L Daubechies, I. Guskov, P. Schroder and 
W. Sweldens, "Wavelets on irregular point sets", Phil. Trans. R. Soc. Lond A., Vol. 357, No. 1760, 
pp. 2397-2413, 1999- For a fixed position X q (m) of the sample points, the irregular subband filtering 

25 operator feat associates D q {m) to C,(m) is a linear operator that we also write O q , Le. D q = O q C q . 

In an exemplary implementation of fee inverse warped subband filtering processor (903), the value of 
D q (m) for m 6 S i>q are computed wife the inverse kernel I m = 0~ T I m , by solving fee system 

Vm€^, £iUOiVO = o. 09) 
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The solution is calculated with a singular value decomposition if the system does not have a unique solutioa 
The inverse subband filtering operator Of 1 for iiregular point sets is then applied on the line values D g (m) 
to recover C q [m). Forallm e S^S^ wesetTTC^L^m)) = C q (m). The insertion system (19) together 
with O" 1 compute the inverse O^ 1 of the linear operator O q previously defined. Any other implementation 
5 of O^ 1 may be used, including a lifting scheme. 

The inverse warped wavelet packet processor (304) in FIG. 3 implements the inverse of the warped 
wavelet packet processor (204) in FIG. 2. It takes in input a waiped wavelet packet transform including 
the warped coefficients at the leaves of a wavelet packet binary tree together with their warping grids, and 
it computes a warped signal. FIG. 10 shows an exemplary configuration of an inverse waiped wavelet 
10 packet processor implemented with Inverse Warped Subband Processors (IWSP) that are distributed along 
the nodes of the same binary tree as the one illustrated in FIG. 5 and used by the waiped wavelet packet 
processor (204). At each node p of this tree, the waiped coefficients WC P and its warping grid WG P 
are reconstructed with an inverse waiped subband processor in a particular direction d, from the children < 
warped coefficients WC 2p and WC^h-i together with their waiping grids WG^ and PTG^+i . The inverse 
15 waiped wavelet packet processor outputs a waiped signal including the waiping grid WG X and the waiped 
coefficients WC\. 

FIG. 9 shows, in block diagram form, an exemplary configuration of an Inverse Waiped Subband Pro- 
cessor (IWSP) along direction d. It implements the inverse of die waiped subband processor illustrated in 
FIG. 8. The processor (901) performs grid union in direction d by taking in input the waiping grids WG^p 
20 and WG^p+i and computing 

WG p (i 7 2 d k) = WChfa k) and WG p (i, 2 d k + l d ) = WG 2p +x(i, k) , 

and RegWG p (t) = 2 d RegWG^i). In (902% the neighborhood connection processor in direction d is 
identical to the neighborhood connection processor (802) in FIG. 8. It computes a waiped neighboihood 
map WN P which is an input together with WG P and die waiped coefficients WC^ and WChm in the 
inverse waiped subband filtering processor (903) along direction d. This processor outputs die waiped 
25 coefficients WC P . Given the waiped neighboihood map WN p and die waiping grid WC P> the processor 
(903) implements die inverse of the waiped subband filtering processor in (803). Its detailed implementation 
was described together with the processor (803). 

Warping Grid Subsa mpling 

The wavelet packet waiping grid subsampling module (302) in FIG. 3 takes in input a waiping grid WG and 
30 computes wavelet packet waiping grids, with a procedure illustrated by the block diagram in FIG. 1 1. For 
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convenience, the input warping grid is relabeled WG\. This warping grid is then subsampled in different 
directions along the same binary tree as the one used to compute the warped wavelet packet processor (204). 
At each node p, the grid subsampling in direction d p is identical to the grid subsampling module in direction 
d that appears in (801) and was previously defined. The wavelet packet warping grids is the family of all 
5 waiping grids WG P at all leaves p 6 P of the binaiy tree. In FIG. 1 1, it corresponds to WG* 9 WG 7 , WGio, 
WGxu WG12 and WGu, i.e. P = {4, 7, 10, 11, 12, 13}. 

Bandeletisation and Inverse Bandeletisation 

The bandeletisation processor (205) takes in input a warped wavelet packet transform composed of warped 
coefficients (WC p ) p eP and corresponding waiping grids (WG p )p^p and outputs bandelet coefficients, by 
10 performing a sequence of 1 -dimensional decorrelations along selected directions on the warped wavelet 
packet coefficients to reduce identified correlations in the signal and obtain a new, sparser representation of 
the signal. Fig. 12 illustrates an exemplary embodiment of a bandeletisation processor. For each input warp- 
ing grid, the direction selection processor (1201) compute a decorrelation descriptor DirWG p (s) depicting 
what land bfaecbifelation islo* be pefformea*dn a region % of waiping grids WG P . The input warping grids 
15 WG P and warped coefficient arrays WC P are split in module (1202) into subgrids WGpa and subarrays 
WCpa wherein the decorrelation descriptor is uniform and equal to some DirWGpo-. The subarrays WCpc 
are then transformed in die bandeletisation along DirWGp^ (1203) into a set of bandelet coefficients. 

The inverse bandeletisation processor (303) takes in input wavelet packet warping grids (WGjOpgP 
and bandelet coefficients and outputs a warped wavelet packet transform composed of computed warped 
20 wavelet packet coefficients {WC p )p£P, and of input wavelet packet waiping grids {WG p )pep. Fig. 13 
illustrates an exemplary embodiment of an inverse bandeletisation processor. Module (1301) computes 
decorrelation descriptors exactly as in module (1201). Module (1302) splits grids WG P into grids WG^ 
exacdy as in module (1202). The inverse bandeletisation along DirWCp* (1303) computes from input 
bandelet coefficients warped coefficient arrays WCpa, which are then merged back into arrays WC P in the 
25 warped grid merging processor (1304). 

In an exemplary embodiment, the regularity descriptor of a grid RegWG p (t) is a vector of n reals 
(bi , . . . , bn) that are dyadic or zero, and indicates what kind of warped coefficient regularity is expected for 
each grid WG P and each region i. For each direction d, bd = 0 means that no regularity has been detected 
along direction d, 0 < 6<f < 1 means that fee warped coefficient array WG P has slow variations along 
30 direction d> and bd > 1 means that the signal has close to periodic variations along direction d, of period 6<f . 
Depending on this regularity descriptor, a decorrelation descriptor DirWG p (i) is computed in module 



27 



WO 2004/056120 



PCT/EP2002/014903 



(1201) for each region i and each wiping grid WG P . The decorrelation descriptor is a vector of n integers, 
(fix , . . . , fa), indicating for each direction d if a decorrelation operation has to be performed along direction 
d (when fa > 0), and then what kind of decorrelation has to be performed along that direction d in module 
(1203). Hence when fa = 0, no decorrelation along direction d is to be performed. When fa « 1, a 
5 decorrelation operator is applied along direction d, like for example a wavelet or wavelet packet transform, 
or a discrete cosine transform. When fa > 1, a decorrelation operator performing a periodic decorrela- 
tion of period fa is applied hi an exemplary embodiment a decorrelation operator performing a periodic 
decorrelation of period fa for a 1 -dimensional array of coefficients (L(m)) m ^ is implemented with a 
one-dimensional wavelet or wavelet packet transform, or a discrete cosine transform applied to subarrays of 
10 L of stride fa, i.e. on the subarrays (L(0 d m + r))a, m+r6 Af for r = 0, . . . , fa - 1. 

In an exemplary embodiment, the decorrelation descriptor DirWG p (i) is computed from RegWG p («) = 
(fa,. -A) by setting DirWG,(t) = {fa,... ,fa) where fa « 0 if b d = 0, fa = 1 if 0 < b d < 1, and 
fa = b d otherwise. In a degenerate case where RegWG p (i) m (0, . . . , 0). the decorrelation descriptor is 
equal to (0, . . . , 0), in which case the bandeletisation performs no decorrelation and outputs coefficients that 
15 are the input-warped waveletpacket coefficients. 

In yet another exemplary embodiment wherein the input warped wavelet packet coefficients are warped 
wavelet coefficients, the decorrelation descriptor DirWGj,(t) is computed in the same way as above, except 
when 0 < b d < 1, in which case fa = 1 if the waiped wavelet coefficients of WC P are low-pass along 
the direction d and fa = 0 otherwise. Warped wavelet coefficients WC P are said to be low-pass along 
20 the direction d if in the chain of warped subband processings to compute WC P from WCu WC P is the 
low-pass output of the last warped subband processing in direction d, or if it has been computed from this 
low-pass output with further warped subband processings in other directions but not in the direction d. 

In yet .another exemplary embodiment, the decorrelation descriptor DirWG p (t) is computed accord- 
ing to one of the two above embodiments, except when all 0 < b d < 1 for all d t in which case we set 
25 DirWGp(i) = (0,-- ,0) instead This corresponds to the case where no further decorrelation is necessary 
because the waiped signal is uniformly regular in all directions and hence the the warped wavelet packet co- 
efficients array is already sparse. The same descriptors are needed in the inverse bandeletisation processor, 
so module (1301) is identical to (1201). 

In an exemplary embodiment, the bandeletisation processor applies decorrelation operators separately 
30 on each region t. In a preferred embodiment, the bandeletisation processor operates across regions when- 
ever this is possible, i.e. whenever two regions are processed according to the same decorrelation descriptor 
DirWGp(i). An exemplary method to apply different decorrelation operators for warping grids correspond- 
ing to different regions i consists in first splitting the input warping grids WG P and warped coefficient 
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arrays WC P into grids (WGp*)^ a*d subarrays {WCpa) a ^ over which the decorrelation descriptor 
DirWGp(t) is the same, and then computing neighborhood connection graphs WN^ for each result- 
ing grid WGp*. In module (1202), WG^ and WC^ are then defined by WG^k) = WG p (i,k) 
if DirWGp(i) = a and WG^fak) = nil otherwise, and similarly WCpe(i,k) = WC p {i>k) if 
5 DirWG p (i) = o and WCpa (i, k) = nil otherwise, where the a index is a decorrelation descriptor. The case 
where the decorrelation operators act separately on each grid is implemented by splitting the warping grids 
WG P and warped coefficient arrays WC P according to the region index. We set WGpa(i, k) = WG p {i, k) 
if a = i and WGpa{i % k) - nil otherwise, and similarly for WCp*, where a is now a region index. In 
all cases, the decorrelation descriptor DirWG p (t) is the same for all region indexes i of the warping grids 
10 WGpa, (i.e. all i such that there exists a k with WG^i, k) ^ nil). We thus write this decorrelation de- 
scriptor DirWG|ar- The decorrelation dimension DimWGp of the wavelet packet warping grids {WG p ) pe p 
is then defined as the maximum number of nonzero entries in a vector DirWGp* for all possible p and <r. 
It is the maximum number of directions along which a bandelet decorrelation will be performed within a 
single grid. 

15 In an inverse bandeletisation, die warping grids WG^ are computed from the warping grids WG P in 
module (1302) in the same way as in module (1202). Module (1302) differs however from module (1202), 
in that it does not compute and output the warping coefficients WC^. 

In the degenerated case where all vectors DirWGp* are zero vectors (0, > . . , 0), then the decorrelation 
dimension DimWGp is equal to 0 and the bandeletisation outputs bandelet coefficients equal to die input 

20 warped wavelet packet coefficients. 

If die decorrelation dimension DimWGp is equal to 1 then the number of directions along which a 
decorrelation operator has to be applied is at most equal to 1, i.e. each vector DirWGpa has at most one 
nonzero entry. This happens most often for 1 -dimensional or 2-dimensional input signals. In this case, 
the decorrelation operator operates along lines of the warped wavelet packet coefficients, and in particular 

25 inside each grid along a single direction. This particular implementation of a bandeletisation is called a 
one-dimensional bandeletisation. For each set of warping grids WG^ the vector DirWGp* either contains 
a single coefficient fa > 0, in which case the corresponding d is denoted d^ and the corresponding value 
0 d is denoted fa, or DirWG^ only contains zero entries, in which case we set fa = fa « nfl. The 
bandeletisation along DirWGj^ in module (1203) then consists in applying when d^ nil a decorrelation 

30 operator along separate lines of each warped coefficient array WCpa along the single direction fa. In an 
exemplary embodiment, the module (1203) computes for each set of warping grids WGpa a neighborhood 
connection graph WNpa(i, fc, ±1, fa) as done for grids WG P in module (802). Then, each set of warping 
grids WGpa is decomposed into a union of connected lines along direction dpa. To do this, we traverse 
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WNpa to build an array of line beginning points Begin^q) = (i, k) that are points ofWG^ that have no 
left neighbor TFJVft fe » = Each line indexed by q of WG^ commences with L m {Q) = 

^ line 58 itera^vety augmented with the right neighbor of the preceding point if (», k) = 
L i*rq( m )> Lp*q{™> + 1) = WJV(«, &, 1, d^), until that point is ml 
5 Then, each line of coefficients C m (ro) = WC^L^m)) for p and a such that ^ = 1 is trans- 
formed into a new array D w with a one-dimensional invertible decorrelation operator. In an exemplary 
embodiment, this invertible decorrelation operator may be a wavelet transform, a wavelet packet transform, 
a discrete cosine or sine transform, computed with a fast algorithm as described in S. Mallat, "A wavelet 
tour of signal processing", Academic Press, 1998. Each line of coefficients for which fi^ > 1 is trans- 
10 formed with a one-dimensional decorrelation operator performing a periodic decorrelation of period ^ va . In 
an exemplary embodiment of a periodic decorrelation operator, for each r = 0, . . . , /fljp* - 1, the subarray 
(Cpvqifipam + r)) m is transformed by the one-dimensional invertible decorrelation operator; and the result 
is stored into {Dpaqifipatn + r)) m . The transformed coefficients of are then stored back into arrays 
WCprt by setting WC^L^m)) = D^ q {m). When = nil, the array WCp^o is equal to WCp*. 
15 -The-bandeleUransform of the. warped coefficient arrays (fVC p ) p ^p is then composed of the coefficient 
arrays WCp^o for all p and a. 

If the decorrelation dimension DimWGp is equal to 1 then the inverse bandeletisation is called an 
inverse one-dimensional bandeletisation and it implements the inverse transform of the one-dimensional 
bandeletisation. In this case, the inverse bandeletisation along DirWGp^ of (1303) is implemented with the 
20 same procedure as (1203), but by replacing die one-dimensional decorrelation operators used in the ban- 
deletisation along DirWGjv by their inverses, called one-dimensional inverse decorrelation operators. In 
the exemplary embodiment, where the decorrelation operators are either wavelet transforms, wavelet packet 
transforms, discrete cosine or sine transforms, their inverse is computed with fast algorithms described in 
S. Mallat, **A wavelet tour of signal processing", Academic Press, 1998. The inverses of one-dimensional 
25 decorrelation operators performing a periodic decorrelation are called periodic inverse decorrelation opera- 
tors. 

We now describe the implementation of a bandeletisation when die decorrelation dimension DimWGp 
is stricdy laiger than 1. In this case there is at least one region where the bandeletisation applies one- 
dimensional decorrelation operators along different directions inside the same region. Inside any such re- 
gion, each decorrelation operator must then be phase-aligned. The bandeletisation is then called a phase- 
aligned bandeletisation. In a preferred embodiment, the bandeletisation along DirWGp* is computed with - 
phase aligned warped wavelet packet transforms, which are applied separately to each pair WGp? and 
WCp<j. Each one-dimensional decorrelation operator is then implemented with a phase-aligned warped 
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, subband processing. For each p and a, we write DirWGp* as OW, . . . jjSp™). To implement periodic 
decollations with phase aligned warped wavelet packet transforms, in an exemplary embodiment the warp- 
ing grids WGpa and waiped coefficients WC^ are subsampled along each direction d for which fipvd > 1. 
To include the non-periodic decorrelation in the same framework, the grids are subsampled by a factor of 
ppad iffiptrd > 1 and not subsampled, or equivalently subsampled by a factor 1, if ft^ = 0 or 1. We thus 
define fi^d = 1 if Ppvd = 0 and jS^d — Ppod otherwise. For each p and a 9 we define Rp^ as die set of 
integer vectors r = (ri,...,r n ) such that 0 < Td<Ppcd ford = l,...,n. The subsampled warping grids 
and warped coefficient arrays WG^ and WCpar are then defined for each r 6 Rpc as 

WG^ii, kx, ...,*„) = WGpa{i, hi A + n, . . Kfa + r„) 
IFC^r (<, fci, . . . , *«) = VTC^Ci, fei A + n, . . . , fenA* + r n ) 

This subsampling for each p, a yields a family of warping grids {WG^peP^e^reR^ and warped coef- 
ficient anays (WC^peP^z^reRp* • Each pair WGjwr* then undergoes a phase warped wavelet 
packet transform exemplified in Fig. 5, with a binary tree T v%tt of one-dimensional warped subband pn>- 
cessirigslo ^ compul^TTG^^^eQ^" ^^{WCpcrqYqiQ^ • Each' binary treeT^includes one-dimensional 
5 decorrelation operators that are phase aligned warped subband processings performed only in directions d 
for which Pp^d > 0. In an exemplary embodiment, die binary tree is a warped wavelet transform binary 
tree over the subset of directions d for which fi^ > 0. When a coefficient array WCpcr is subsampled 
along direction d with Bfxrd > 0, die phase aligned warped subband processing along direction d is said 
to perform a periodic decorrelation. The bandelet coefficients, which are die output of the bandeletisation 
10 illustrated in Fig. 12 are the coefficients (WC pa rq)p€P 1 ae^ t reR J H f1 q€Qpa» These are die input of the inverse 
bandeletisation illustrated in Fig. 13. 

The inverse phase-aligned bandeletisation implements the inverse of the operator implemented by die 
phase-aligned bandeletisation processor; hi this case, the inverse bandeletisation along DirWGpo- (1303) 
computes from the waiping grids WGpe all the subsampled warping grids WGp^r with die same procedure 
as in die module (1203). Then for each p, cr 7 r, the warping grids {WCparqiqeQ^ are computed from 
WGpar with a tree of grid subsamplings corresponding to T^, as illustrated in Fig. 1 1. For each p, <r, r, the 
warping grids (WG P arq)qeQ pir and warped coefficient arrays (WCp ar q)qeQ pir are transformed by a phase 
aligned inverse warped wavelet packet processor implemented with die binary tree Tpa of phase aligned 
inverse warped subband processings, to compute a waiped coefficient array WCpcr. When Pp^d > 0, 
the phase aligned inverse warped subband processings are said to perform a periodic inverse decorrelation. 
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Then, for each p and a, the warped coefficient arrays {WC^reR^ are merged back into WCp* by setting 

WCprii, kt ft + r lf . . M kn$n + r n ) = WCprrih *1, • • - , *n) 

for each r in and fc such that WG^kiPi + ri,.-.,fc„& + r n ) ^niL 

For both the inverse one-dimensional bandeletisation processor or die inverse phase-aligned bandeleti- 
sation processor, for each p € P the waiped coefficient arrays {WCpa) ae v are merged back into WC P 
in the warping grid merging processor (1304) by setting WC p (i,k) = WC^i, k) if a is such that 

It is apparent to those skilled in the art that is also possible to combine one-dimensiona] bandeletisa- 
tion and phase-aligned bandeletisation. Denoting DimWG^ the decollation dimension of WGpv, which 
is the number of nonzero entries in DirWGp*, the combined bandeletisation consists in applying a one- 
dimensional bandeletisation on all grids WG^ such that DimWGp, = 1, and applying a phase-aligned 
10 bandeletisation on all grids WG^ such that DimWGj*, > 1. The same rule is then used for the inverse 
combined bandeletisation. 

Restoration System 

A restoration system is implemented by using a restoration processor in die processor (102) of FIG* 1. 
A restoration processor computes a processed warping geometry and processed bandelet coefficients, and 

15 an inverse warped wavelet packet bandelet processor takes in input the processed warping geometry and 
processed bandelet coefficients to reconstruct a restored signal. 

In an exemplary embodiment, a restoration processor removes additive noises, by applying diagonal 
operators to bandelet coefficients. A diagonal operator applied to a bandelet coefficients stored in WC q (i , k) 
computes 6 q (WC q (i, k)) where 6 q (x) is a linear or non-linear function. The function 9 q can be chosen to 

20 be a thresholding function which sets to zero the bandelet coefficients whose amplitude are smaller than a 
threshold T q whose value adapted to die noise properties. Examples of state-of-the-art linear or thresholding 
estimators are described in S. Mallat, A Wavelet Tour of Signal Processing, 2nd edition. Academic Press, 
San Diego, 1999. A rpgularization operator can also be applied to the warping geometry calculated from the 
input noisy signal, in order to suppress the effect of die noise on the signal geometry. 

25 A restoration processor can also implement a signal deconvolution system, for example to debluor 2- 
dimensionaJ images. In an exemplary embodiment, a diagonal operator is applied to bandelet coefficients 
where the diagonal operator is designed to approximate the deconvolution filter by multiplying each ban- 
delet coefficient with an appropriate constant A restoration processor can also implement any state of die 
operators on the bandelet coefficients and the warping geometry to restore an output signal for restoration 
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systems, including inverse scattering systems, tomographic reconstruction systems, interpolation systems, 
or supeiresolution systems which reconstruct higher resolution signals. 

Compression System 

FIG. 14 illustrates an exemplary implementation of a signal compression system using the warped wavelet 
5 packet bandelet processor illustrated in FIG. 2. The input signal is transformed in (1401) into a warping 
geometry and a set of bandelet coefficients. The warping geometry is computed by a geometrical segmenta- 
tion module (201) in FIG. 2, which includes a quantizer to compute the deformation parameters that specify 
the geometrical deformations in each signal region. The quantizer (1402) outputs quantized bandelet coeffi- 
cients that are encoded together with the warping geometry by a binary coder (1403) to produce a bitstream. 

10 The quantizer (1402) may be any state of the art scalar or vector quantizer, adapted to the type of input signal 
and the configuration of the warped wavelet packet bandelet processor. In an exemplary embodiment, fee 
quantizer is a uniform quantizer of a constant bin size T and with a zero bin twice larger than the other bins, 
and equal to [-T,T]. The binary coder (1403) may be any state of fee art binary coder, including entropy 
coders such as an arithmetic coder, or a Huffman coder, or any entropy coder, as in T. M. Cover and J. A. 

15 Thomas, Elements of Information Theory, Wiley Series in Telecommunications, John Wiley & Sons, 1991. 
The entropy coding may also use multiple contexts to further reduce fee code size. 

An exemplary embodiment of such a compression system when n = 1 is a compression system for 
electro-cardiograms (ECG), where fee periodic decollation in the bandeletisation is particularly pertinent 
for deconrelating the periodic structure of the ECG. Yet another exemplary embodiment of such a compres- 

20 sion system when n = 2 is an image compression system computed wife a warped wavelet transform, wife 
a warping geometry feat defines regions that are unions of rectangles. For particular class of images such 
as human faces, the warping geometry and its coding takes advantage of prior information known on this 
class of images. In yet another exemplary embodiment for special classes of images such as fingerprint 
images, fee warped wavelet packet transform is adapted to fee property of this class of images. The periodic 
25 decorrelation in fee bandeletisation is pertinent for decorrelating fee periodic ridge structures existing in 
fingerprint images. 

Yet another embodiment of such a compression system, when n = 3 is a video compression system 
using a warping geometry that includes a directional geometric deformation obtained by applying time 
displacements to 2-dimensional geometrical deformations using estimated motion vectors. Yet another ex* 
30 emplary embodiment when n — 3 of such a compression system is a compressor for seismic 3-dimensional 
volumetric data or for 3-dimensional medical data. Yet another embodiment of such a compression system 
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when n = 4 is a compressor for light field volumes of data, such as the data obtained by composing digital 
pictures of an object taken with varying tilt and pan angles. 

FIG. 15 illustrates an exemplary implementation of a signal decompression system using the inverse 
warped wavelet packet and bandelet transform illustrated in FIG. 3. The processor (1501) decodes the input 
5 bitstream and outputs the warping geometry and the quantized bandelet coefficients, from which the inverse 
warped wavelet packet bandelet processor (1502), illustrated in FIG. 3, computes an output reconstructed 
signal. 

Feature Extraction System 

The present invention includes a system that computes a signal feature vector from the warping geometry 
10 and the bandelet coefficients, for pattern recognition applications including content based signal indexing 
and retrieval from a database, signal matching, or for detection and classification. The input signal is trans- 
formed into a warping geometry and bandelet coefficients by the warped wavelet packet bandelet processor 
illustrated in FIG. 2. 

In a first configuration of the system illustrated in FIG. 16, the warping geometry and the bandelet co- 
15 efficients of an input signal are computed with a warped wavelet packet bandelet processor (1601) and a 
feature vector is calculated from the warping geometry and the bandelet coefficients by using a state of 
the art feature extraction procedure (1602). In an exemplary embodiment, the feature vector is used for 
a content-based indexing system. The feature extraction can be computed with histogram techniques ap- 
plied to transformed bandelet coefficients as in M. K. Mandal and T. Aboulnasr, "Fast wavelet histogram 

20 techniques for image indexing", Computer Vision and Image Understanding, vol 75, no. 1/2, pp. 99-110. 
Bandelet coefficients may be equal to warped wavelet coefficients for a degenerated bandeletisation where 
the decorrelation dimension is 0. The feature parameters obtained from the warping geometry describe geo- 
metrical signal properties whereas the parameters obtained from bandelet coefficients describe the evolution 
of the signal values along this geometry. For a content based signal retrieval from a data basis of signals, 

25 the signal feature vector is compared with the feature vectors of the signals stored in a data basis, using 
an appropriate distance measure. In yet another exemplary embodiment, feature vectors are used for signal 
matching and recognition by comparing the respective feature vectors of two signals. The feature extraction 
is performed with any state of the art procedure adapted to matching applications for a particular class of 
signals, for example fingerprint images or human face images. In yet another exemplary embodiment, fea- 

30 ture vectors are used to classify signals. For an electro-cardiogram, a feature vector can be used by a state 
of the art classifier to make a diagnosis on the electro-cardiogram signal and detect anomalies. 



34 



WO 2004/056120 



PO7EP2002/014903 



To reduce bandwidth and storage requirements, signals are represented in a compressed form. In a 
second exemplary configuration of the system illustrated in FIG. 17, the signal feature vector is computed 
(1701) from the compressed bitstream produced by the binary coder (1403) of the warped wavelet packet 
bandelet compression system illustrated in FIG. 14. Any state of the art technique may be used to compute 
the feature extraction. For content-based indexing and retrieval from a database, in an exemplary embodi- 
ment, the parameters of the feature extraction are obtained with state of the art histogram procedures as in E. 
Feig and C.-S. Li, Computing image histogram from compressed data," Proceedings of SPEE 2898, 1996, 
pp. 1 18-124. For matching or classification of signals such as fingerprint images, images of human feces or 
electro-cardiograms or any other type of signals, any state of the art feature extraction system can be used 
to compute the feature vector: 

While a detailed description of presently exemplary embodiments of the invention has been given above, 
various alternatives, modifications, and equivalents will be apparent to those skilled in the art For example, 
while different components of the warped wavelet packet processor and inverse warped wavelet packet 
processor of the present invention are described herein as performing certain specific functions, one skilled 
in the art will appreciate that other components or circuits in the service module may perform some or all of 
such functions without varying from the spirit of the invention. Therefore, the above description should not 
be taken as limiting the scope of the invention which is defined by the appended claims. 
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